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Abstract. For chaotic cavities with scattering leads attached, transport properties 
can be approximated in terms of the classical trajectories which enter and exit the 
system. With a semiclassical treatment involving fine correlations between such 
trajectories we develop a diagrammatic technique to calculate the moments of various 
transport quantities. Namely, we find the moments of the transmission and reflection 
eigenvalues for systems with and without time reversal symmetry. We also derive 
related quantities involving an energy dependence: the moments of the Wigner delay 
times and the density of states of chaotic Andreev billiards, where we find that the gap 
in the density persists when subleading corrections are included. Finally, we show how 
to adapt our techniques to non-linear statistics by calculating the correlation between 
transport moments. 

In each setting, the answer for the n-th moment is obtained for arbitrary n (in 
the form of a moment generating function) and for up to the three leading orders in 
terms of the inverse channel number. Our results suggest patterns which should hold 
for further corrections and by matching with the low order moments available from 
random matrix theory we derive likely higher order generating functions. 
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1. Introduction 

Transport through a chaotic cavity is usually studied through a scattering description. 
For a chaotic cavity attached to two leads with Ni and N2 channels respectively, the 
scattering matrix is an iV x iV unitary matrix, where = A'^i + It can be separated 
into transmission and reflection subblocks 

which encode the dynamics of the system and the relation between the incoming 
and outgoing wavefunctions in the leads. The unitarity of the scattering matrix 
S'^S = I = 55^ leads to the following relations (among others) 

rln + tH = In, , r2rl + tt^ = 1^, , (2) 

while the transport statistics themselves are related to the terms in (2) involving the 
scattering matrix (or its transmitting and reflecting subblocks) and their transpose 
conjugate. For example, the conductance is proportional to the trace Tr ^tH'j (Landauer- 
Biittiker formula [1, 2, 3]), while other physical properties are expressible through higher 
moments like Tr 

There are two main approaches to studying the transport statistics in clean ballistic 
systems: a random matrix theory (RMT) approach, which argues that S can be 
viewed as a random matrix from a suitable ensemble, and a semiclassical approach 
that approximates elements of the matrix S by sums over open scattering trajectories 
through the cavity. 

It was shown by Bliimel and Smilansky [4, 5] that the scattering matrix of a 
chaotic cavity is well modelled by the Dyson's circular ensemble of random matrices 
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of suitable symmetry. Thus, transport properties of chaotic cavities are often treated 
by replacing the scattering matrix with a random one (see [6] for a review). The 
eigenvalues of the transmission matrix tH then follow a joint probability distribution, 
which depends on whether the system has time reversal symmetry or not, and from 
which transport moments and other quantities can be derived. Though the conductance 
and its variance were known for arbitrary channel number [7, 8], other quantities were 
limited to a diagrammatic expansions in inverse channel number, see [9]. However, the 
RMT treatment has recently experienced a resurgence due to the connection to the 
Selberg integral noticed in [10]. Following the semiclassical result for the shot noise [11], 
the authors of [10] used recursion relations derived from the Selberg integral to calculate 
the shot noise and then later all the various moments up to fourth order for arbitrary 
channel number [12]. 

Since then a range of transport quantities have been treated, for example the 
moments of the transmission eigenvalues for chaotic systems without time reversal 
symmetry (the unitary random matrix ensemble) [13, 14, 15, 16] and those with time 
reversal symmetry (the orthogonal random matrix ensemble) [15, 16]. For the unitary 
ensemble, the moments of the conductance itself were also obtained in [14] and, using a 
different approach, in [17] which was later extended to the moments of the shot noise [18]. 
Building again on the Selberg integral approach, the moments of the conductance and 
shot noise have been derived for both symmetry classes [19]. Interestingly these results, 
though all exact for arbitrary channel number, are given by different combinatorial sums, 
and the question of how they are related to each other is still open in many cases. 

On the other hand, the semiclassical approach makes use of the following 
approximation for the scattering matrix elements [20, 21, 22] 

^ Yl ME)e'^'^^^\ (3) 

which involves the open trajectories 7 which start in channel i and end in channel o, 
with their action 5*^ and stability amplitude A^. The prefactor also involves which 
is the average dwell time, or time the trajectory spends inside the cavity. For transport 
moments we consider quantities of the type 

<^'[^"^']">~((^n^ ^ ^-„^;;e^''--";'\ . (4) 

where the trace means we identify in+i = ii and where X is either the transmitting 
or the reflecting subblock of the scattering matrix. The averaging is performed over a 
window of energies E. 

The choice of the subblock X affects the sums over the possible incoming 
and outgoing channels, but not the trajectory structure which involves 2n classical 
trajectories connecting channels. Of these, n trajectories 7^, j = l,...,n, contribute 
with positive action while n trajectories 7^- contribute with negative action. In the 
semiclassical limit of ^ — )■ we require that these sums cancel on the scale of h so that 



Figure 1. (a) The semiclassical trajectories for the second moment travel around a 
closed cycle. By collapsing the trajectories onto each other as in (b) we can create 
a small action difference and a trajectory quadruplet with a single encounter that 
contributes to the shot noise at leading order in inverse channel number, (c) For the 
conductance, a trajectory pair with a similar encounter provides the first subleading 
order correction for systems with time reversal symmetry. 

the corresponding trajectories can contribute consistently when we apply the averaging 
in (4). 

The main idea of the semiclassical treatment is that, in order to achieve a small 
action difference, the trajectories {7^}, must follow the path of trajectories {7-,} most 
of the time, deviating only in small regions called encounters. This is best illustrated 
with an example. In figure 1(a) the schematic depiction of the trajectories is shown 
for the case n = 2. We have 2 trajectories 7 depicted by solid lines, 71 : ii — )■ Oi 
and 72 : ^2 — ^ 02, and 2 trajectories 7' depicted by dashed lines, 7^ : «2 — ^ Oi and 
72 : — )■ 02. Figure 1(b) shows one possible configuration for achieving a small action 
difference: trajectory 7^ departs from the incoming channel 12 following the path of 
trajectory 72. Then, when the trajectories 71 and 72 come close to each other in phase 
space (thus the term 'encounter'), trajectory 7^ switches from following 72 to following 
7i, before arriving at its destination channel oi. The trajectory 72 does the opposite. 
The picture in figure 1(b) is referred to as a 'diagram'; it describes the topological 
configuration of the trajectories in question, while leaving out metric details. The task 
of semiclassical evaluation can therefore be divided into two parts: evaluation of the 
contribution of a given diagram by integrating over all possible trajectories of given 
structure and enumeration of all possible diagrams. 

Historically, the semiclassical treatment started with the mean conductance 
(Tr [t^t]), involving a single trajectory and its partner. The leading order contribution 
comes from trajectory pairs that are identical — the so-called diagonal approximation 
which was evaluated in [23, 24]. The first non-diagonal pair was treated in [22] and 
involved a single encounter where one trajectory had a self-crossing while the partner 
avoided crossing as in figure 1(c). Such a pair can only exist when the system has time 
reversal symmetry and its contribution was shown to be one order higher in inverse 
channel number, 1/N, than the diagonal terms. The expansion to all orders in inverse 
channel number was then performed in [25] for systems with and without time reversal 
symmetry by considering arbitrarily many encounters each involving arbitrarily many 
trajectory stretches. 

Importantly, the work on the conductance [25] showed that the semiclassical 
contribution of a diagram can be decomposed into a product over the constituent parts of 
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the diagram, greatly simplifying the resulting sums. In fact for the second moment, the 
shot noise, all such diagrams were generated in [11] and with them the full expansion in 
inverse channel number. This, along with the conductance variance and other transport 
correlation functions, as well as the semiclassical background was covered in detail in 
[26]. 

However, the method for diagram enumeration considered in [25, 11, 26] becomes 
unwieldy for higher moments, which encode finer transport statistics. To the leading 
order in 1/A^, the higher moments were derived in [27]. The semiclassical approach 
requires a large number of channels in each lead, 1 <^ Ni,N2, but to unambiguously 
separate the orders in inverse channel number one may additionally assume that both 
A^^i and N2 are of the same order as N. For example, the result of [27] was in terms 
of the variable ^ = NiN2/N'^ which should then be constant and introduce no further 
channel number scaling. We therefore make the same assumption in this article when 
describing the different orders in 1/A^, though of course a different scaling, say keeping 
Ni fixed so that ^ ~ 0(1/A^), may simply lead to a mixing of the different 'orders' 
without changing the individual results. 

The diagrams contributing at the leading order to the n-th moment were shown 
in [27] to be trees. The tree expansions turned out to be very well suited to analysis 
of other interesting physical quantities, such as the statistics of the Wigner delay times 
[28], which are a measure of the time spent in the scattering region, and the density 
of states in Andreev billiards [29, 30]. If we imagine replacing the scattering leads by 
a superconductor we have a closed system called an Andreev billiard. Each time an 
electron inside the system hits the superconductor it is reflected as a hole retracing its 
path until it hits the superconductor and is retrorefiected as an electron again. Wave 
interference between these paths leads to significant effects, most notably a complete 
suppression of the density of states for a range of energies around the Fermi energy. 
Similarly, strong effects on the conductance (of the order of the mean conductance) can 
also be seen if we attach additional superconducting leads to our original chaotic cavity 
(making a so-called Andreev dot) [31, 32]. The size of these effects make such systems 
especially interesting for a semiclassical treatment. But treating these effects effectively 
requires knowledge of all the higher moments, and gives us strong reason to go beyond 
low n. 

One particular nicety of the semiclassical approach is that it can incorporate, in 
a natural way, the effect of the Ehrenfest time. This is the time scale that governs 
the transition from classical dynamics to wave interference, which dominates when the 
Ehrenfest time is small (on the scale of the typical dwell time). For larger Ehrenfest 
times, the competition between the different types of behaviour leads to quite striking 
features, like an additional gap both in the density of states of Andreev billiards and in 
the probability distribution of the Wigner delay times [30, 33]. Semiclassically we can 
explicitly track the effect of the Ehrenfest time all the way to the 'classical' limit, which 
can only be achieved using RMT by postulating the Ehrenfest time dependence of the 
scattering matrix. 
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Alongside the case of ballistic systems, the typical chaotic behaviour and transport 
statistics can also be induced by introducing disorder in the system. For weak disorder 
the transport properties coincide with those obtained from RMT, and one can also obtain 
the full counting statistics at leading order, as well as weak localisation corrections and 
universal conductance fluctuations, using circuit theory [34]. If the disorder averaging 
is treated using diagrammatic perturbation theory (see, for example, [35]) multiple 
scattering events can be summed, in the limit of weak disorder, as a ladder diagram 
known as a Diffuson. This corresponds to the parts of semiclassical diagrams where 
the trajectories are nearly identical, as on the left of the encounter in figure 1(c). The 
disordered systems' counterpart of the loop on the right of the encounter in figure 1(c), 
traversed by trajectories in opposite directions, is called a Cooperon, while the encounter 
itself corresponds to a Hikami box [36]. Although transport properties like the weak 
localisation diagram related to figure 1(c) and conductance and energy level fluctuations 
[37] can be treated diagrammatically, usually powerful field-theoretic methods involving 
the nonlinear a model are used (see [38] for an introduction). These methods can 
treat both weak and stronger disorder non-perturbatively, and by using supersymmetry 
[39, 40] a large range of transport and spectral properties can be obtained, for open 
and closed systems correspondingly. More importantly, the applicability of RMT for 
weakly disordered systems can be justified and RMT shown to be the zero- dimensional 
variant of the a model [41, 42]. Alongside the supersymmetric a model, there is also 
the replica a model which is particularly useful for perturbative expansions. This leads 
to a diagrammatic expansion, with diagrams that can be reinterpreted as correlated 
semiclassical trajectories [43]. In fact this connection between semiclassical diagrams 
and disorder diagrams from the replica a model lay behind the semiclassical treatment 
of energy level correlations in closed systems [44, 45, 46] which in turn led to the 
semiclassical treatment of transport [22, 25, 26] discussed above. 

To summarize, there are established semiclassical tools for the analysis of (4) for 
small n to all orders of 1/A^ and for all n but only to the leading order of 1/A^. It 
is the purpose of this article to start closing this gap. For all n we derive the next 
two corrections for (4) and related quantities. We show that the contributing diagrams 
can be generated by grafting trees onto the 'base diagrams', which can be obtained 
by 'cleaning' the diagrams used in [25]. We therefore first review the leading order 
tree recursions in section 2 before treating transport moments beyond the leading order 
in section 3. We start by cleaning the diagram of figure 1(c) which gives the first 
subleading order orthogonal correction. Grafting trees onto the base diagram leads to a 
generating function, which we apply to calculate the moments of the transmission and 
reflection eigenvalues. Proceeding to the next order in 1/A^ we then treat the second 
subleading order diagrams for the unitary and the orthogonal case. For the moments of 
the reflection and transmission eigenvalues we find that our generating functions simplify 
and become rather straightforward. 

The graphical recursions we use provide a new insight into the leading order terms 
which is particularly useful for energy dependent correlation functions. Such correlation 
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functions are needed for a treatment of the density of states of Andreev billiards, which 
we consider in section 4 where we find that the hard gap, previously found at leading 
order in persists at least for the next two orders. Also derivable from energy 

dependent correlation functions are the moments of the Wigner delay times, treated in 
section 5, and we find that the corrections at each order in are also generated by 
relatively simple functions. Of course, the transport moments in (4) are only one type 
of transport quantity, and we finally look at non-linear statistics in section 6 and see 
how their treatment follows naturally from the previous semiclassical considerations. 

We shall be comparing our semiclassical results with the prediction of RMT, where 
those predictions are available: previously (of the quantities treated here) only the 
moments of the transmission amplitudes for systems without time reversal symmetry 
have been given for an arbitrary number of channels [13, 14]. Explicit results for systems 
with time reversal symmetry have just been derived [15, 16] and we were pleased that 
[15] shared those results with us beforehand. The moments of the Wigner delay times 
for both symmetry classes have also been obtained [15]. 

Of the recent RMT results, it is those concerned with the asymptotic expansion as 
the number of channels increases, currently to leading order [47, 48, 49], that particularly 
connect with the work here. Semiclassically, without the equivalent of the Selberg 
integral, we are still restricted to an expansion in inverse powers of the channel number, 
but as we shall see the semiclassical treatment leads to explicit and surprisingly simple 
generating functions at each order in inverse channel number. This simplicity until now 
remained hidden in the combinatorial sums of the RMT results and may suggest ways 
of simplifying those results and of highlighting the underlying combinatorial structure. 

2. Subtrees 

The semiclassical treatment of the conductance beyond the diagonal contribution, 
starting [22] with the trajectory pair depicted in figure 1(c), required two main 
ingredients. The first was to estimate how often a trajectory would come close to 
itself and have a self-encounter. This is performed using the global ergodicity of the 
chaotic dynamics. The second was that, given such an encounter, we can use the local 
hyperbolicity of the motion to find the partner trajectory which reconnects the stretches 
of the original trajectory in a different way. Then one can determine the action difference 
between the two trajectories and hence their contribution in the semiclassical limit. 
When treating diagrams with more numerous and more complicated encounters, the 
authors of [25] showed that these two ingredients allowed them to express the total 
contribution as a product of integrals over the encounters and over the 'links', the 
trajectory stretches which connect the encounters together. Performing these integrals 
then led to simple rules for the contributions of the constituent parts of any diagram, 
and essentially reduced the problem down to the combinatorial one of finding all the 
possible diagrams. For the first two transport moments, this was done [26] by cutting 
open the periodic orbit pairs that contribute to spectral statistics [44, 45, 46]. 
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For the higher moments, as shown in [27], the diagrams that contribute at leading 
order in inverse channel number are rooted plane trees. The reason is simple: according 
to the semiclassical evaluation rules of [26], every encounter contributes a factor of — 
while every link contributes a factor of 1/N . The leading order is thus achieved by a 
diagram with the minimal possible difference between the number of links (edges) and 
encounters (internal vertices). It is a basic fact of graph theory that this difference 
is minimized by trees; each independent cycle in a graph adds one to this difference. 
Thus to go beyond the leading order one needs to consider diagrams with an increasing 
number of cycles. We will approach this task by describing the topology of the cycles 
using 'base diagrams' — graphs with no vertices of degree 1 or 2 — and then grafting 
subtrees onto the base diagrams. 

Adding a subtree does not change the order of the contribution in inverse channel 
number 1/N but adds more incoming and outgoing channels thus changing the order 
of the moment n. Because we will be joining the trees to existing structures, unlike the 
treatment in [27, 28, 29, 30], here we do not root our trees in an incoming channel, but 
at an arbitrary point. These trees then correspond to the restricted trees in [27, 30] and 
will be referred to as 'subtrees'. We also note that the generating function variables we 
use here have slightly different definitions than in [27, 28, 29, 30]. Our present choice 
is more appropriate for the subleading orders and the different transport quantities 
considered. 

We now summarize the derivation of the subtree generating functions which were 
introduced in [27] and further developed in [28, 30]. A subtree consists of a root, 
several vertices of even degree (called 'nodes', they correspond to encounters between 
various trajectories) and 2n — 1 vertices of degree one (called 'leaves', they correspond 
to incoming or outgoing channels). The leaves are labelled z or o alternatingly as we 
go around the tree anti-clockwise. There are two types of subtrees: the /-subtrees have 
leaves labelled o^, ik+i, o^+i, ik+2 etc. The label ik would correspond to the root if 
we were to label it too. The /-subtrees have leaf labels ik+i-, o^+i, ifc+2, 0^+2 etc. The 
reference index k depends on the location of the subtree on the diagram. 

It is possible that an encounter happens immediately as several trajectories enter 
the cavity from the lead or exit the cavity into the lead. To keep account of these 
situations, we say that an Z-encounter (node of degree 21) may 'i-touch' the lead if it is 
connected directly to / incoming channels (leaves with label i) and 'o-touch' if connected 
to / outgoing channels. When an encounter touches the lead, the edges connecting it 
to the lead get cut off and all the channels must coincide, although in the diagrams we 
keep short 'stubs' to avoid changing the degree of the encounter vertex. 

We define the generating functions f{x,Zi,Zo) and / = f{x,Zo,Zi) which are 
counting /- and /-trees correspondingly. The meaning of the variables x = (x2, X3, . . .), 
Zi = {zi^2, 2i,3, • • •) and Zo = {zo,2, ^0,3, • • •) is as follows: 

• Xi enumerate the /-encounters that do not touch the lead, 

• Zi^i enumerate the /-encounters that z-touch the lead. 
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(a) 



(b) 
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Figure 2. The subtree shown in (a) is cut at its top node (of degree 4) creating 
subtrees (b)-(d). Subtree (c) has the incoming and outgoing directions reversed. The 
lower vertex in (a), and hence (c), is i-touching the lead so that the channels and 
i5 (not shown) coincide (14 = 15). This is represented by the short stubs, and the 
encounter now starts in the incoming lead. 



• Zo,i enumerate the /-encounters that o-touch the lead. 

For example, the coefficient of X3X2Zi^2 gives the number of trees with one 3-encounter 
(a vertex of degree 6) and three 2-encounters (vertices of degree 4), one of which i- 
touches the lead. An example of such a tree is given in figure 2(a). We note that if 
an encounter may touch the lead, the generating function includes (and sums) both 
possibilities: touching and non-touching. For example, the left-most vertex of the tree 
in figure 2(a) may o-touch the lead, but this possibility is counted separately. 

In addition we will use several secondary parameters that will allow us to adapt 
the subtree generating functions to each of the four quantities considered in the paper. 
These parameters are: 

• y is the semiclassical contribution of an edge (link) 

• c,- and Cq cirG the contributions of an incoming and an outgoing channel 

• 0" is a special correction parameter for the situation when an o-touching node is 
directly connected to an i-channel (a = everywhere apart from section 5). 

We obtain a recursion for the functions / and / by cutting the subtree at the top 
encounter node. If this node is of degree 21, this leads to 2/ — 1 further subtrees as 
illustrated in figure 2. Assuming we started with an /-subtree, / of the new subtrees 
also have type /, while the remaining / — 1 are /-subtrees. Thus an /-subtree with an 
/-encounter at the top contributes yxif y^ to the generating function /. Additionally, 
we consider the possibility for the top node of an /-subtree to o-touch. In this case 
its odd-numbered further subtrees are empty stubs and the even-numbered subtrees 
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are still arbitrary, leading to the contribution yZo^i{f + cr)'~^. Here we have included a 
correction term a which is used in section 5 to control the contribution of any /-subtree 
that consists of one edge and directly connects an incoming and outgoing channel and 
is set to in the rest of the paper. 

We start our recursion relation at the value for an empty tree, which consists of a 
link (with the factor y) and an outgoing channel (providing a factor Co), 

oo 

f = yco + yY, [xiff'-' + z^^iif + aY^'] . (5) 

1=2 

The recursion is similar for /, with the roles of i- and o-variables switched, 

oo 

f = yc^ + yY.\xlJ'f-' + z.,lf'']■ (6) 



1=2 

2.1. Reflection 

For the reflection into lead 1 we will consider the generating function 

oo 

^(^) = E^"(Tr[^Kr), (7) 

n=l 

where the power of s counts the order of the moments. For the individual semiclassical 
diagrams we make use of the diagrammatic rules of [26], where each link contributes a 
factor of while each encounter provides the factor — A^. Each channel is in lead 1, 
so can be chosen from the A'^i available and provides this factor. When an encounter 
starts (or ends) in the lead, all the incoming (or outgoing) channels must then coincide 
in the same channel, leading again to the factor A'^i. Bearing in mind the meaning of 
the variables introduced above, we therefore have to make the following semiclassical 
substitutions: 

y = ^, xi = -N, Zi^i = Zo,i = r^Ni, Ci = Co = rNi, a = 0, (8) 

where we have introduced r whose power counts the total number of channels and 
which allows us to keep track of the total contributions to different moments. The n-th 
moment involves 2n channels so we have the relation s = r^. Each channel factor c then 
includes the factor r, while the formula for Zi^i in (8) accounts for the fact that when 
an /-encounter enters the incoming channels we have / channels coinciding but only a 
single channel factor. 

If we define (i = Ni/N, the subtree recursions (5) and (6) both become 

oo oo 

/ = <i-E/''"'+<iE^'"V'-^ (9) 

1=2 1=2 

Performing the sums (where the terms / and r(i correspond to / = 1 of the sums) this 
is 

= h -!-^L_ (10) 
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which can be written as the quadratic 



1 - ^1 - 4^r2 



r(l-Ci)/'-/ + rCi = 0, f= ^ 2r(l-CiT ' ^^^^ 
where ^ = — Ci) ^-^id where we take the solution whose expansion agrees with the 
contributions of the semiclassical diagrams. 

2.2. Transmission 

For the transmission we treat the function 

oo 

T(.) = 5^."(Tr[ttt]"), (12) 

n=l 

and to distinguish it more clearly from the reflection we will call the corresponding 
subtree generating function / = here. For the transmission, the equations are a bit 
more complicated than for the reflection because 7^ in general. For the substitution 
we need 

y = ^, xi = -N, Zi^i = r^Ni, Zo,i = r'iVa, 

c^ = rNi, Co = rN2, a = 0, (13) 

where the only difference from (8) is that the outgoing channels are now in lead 2 and 
can be chosen from the N2 available. 

The contribution of the subtrees (5) once summed becomes 

0= ^ + ^^, or = rC2 + rCi00, (14) 

1 — 00 1 — r0 

with (2 = N2/N and using that Ci + C2 = 1- Likewise (6) becomes 

r(^i 

= -J + -, or = rCi + rC200, (15) 

1-/0 1 - r0 

where, as before ^ = ^i(l — ^1) = (1(2- With /i = 00 we have 

r^^h^ + [r\l-2O-l]h + r^^ = 0, (16) 

from which we can find the equations 

0'- (^^^^ + 2r^ + 1 = 0, 0'- (^i^^ + 2r^ + 1 = 0. (17) 



3. Transport moments 

By a simple counting argument, the order of a diagram in terms of inverse channel 
number is the number of edges minus the number of vertices (both leaves and nodes). 
Thus a diagram contributes at the order {1/N)^~^, where /3 is the number of the 
independent cycles in the diagram (also known as the cyclomatic number or the first 
Betti number, hence the notation). The leading contribution thus comes from tree 
diagrams which have /3 = and the next contribution comes from diagrams with one 
cycle. 
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(a) ii 




(b) 




«2 



Figure 3. The correlated trajectory quadruplet in (a) which contributes to the second 
moment at leading order in inverse channel number can be redrawn as the ribbon tree 
in (b) by 'untwisting' the encounter. The four trajectories themselves can be read off 
from the boundary walk shown. At subleading order in inverse channel number, we 
start with the correlated periodic orbit pair in (c) which can be represented as the 
graph in (d) with corresponding boundary walks. Cutting the periodic orbit along the 
left link (which is traversed in the same direction by the orbit and its partner) creates 
the correlated trajectory pair in (e) which contributes to the first moment. Changing 
this diagram into a graph we arrive at the structure in (f ) which is a Mobius strip with 
an empty subtree inside and outside the loop. The intertwined S's in diagrams (d) and 
(f) represent twists in the corresponding ribbon links. 



3.1. First orthogonal correction 

A diagram with one cycle can be thought of as a loop with trees grafted on it. But there 
is a twist. The reconstruction of the trajectories' structure from a tree, see [27], was 
done by means of the boundary walk. It helps to visualise the edges of the tree as strips, 
a model that is called a fat or ribbon graph in combinatorics. This fixes the circular 
order of edges around each vertex and, going along the boundary, prescribes a unique 
way to continue the walk around a vertex (see [50] for an accessible introduction). The 
trajectories 7j of equation (4) are then read as the portions of the walk going from ij 
to Oj. The trajectories 7^, on the other hand, appear in reverse as portions of the walk 
going from Oj to ij+i- For example, the diagram in figure 3(a) which contributes at 
leading order in 1/iV can be redrawn as the tree in figure 3(b) with the corresponding 
boundary walk shown. 

The trace in (4) means that the boundary walk is closed and the equality of total 
actions implies that each edge of the diagram is traversed twice (once by 7 and once 
by 7'). This means that a valid diagram must have one face. In particular, there must 
be a way for the walk to cross from inside to outside of the cycle of the first correction 
diagram. The diagram thus has the topology of a Mobius strip with (ribbon) trees 
grafted on the edges. We will refer to the diagram without any trees (Mobius strip in 
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Figure 4. To obtain the base structure in (a) we can simply remove the empty subtrees 
of figure 3(f). Appending subtrees to (a) we can then create all the possible graphs, 
but for the graph to remain a Mobius strip we need an odd number of odd nodes, as for 
example in the graph in (b). We draw the boundary walk in (c) where we truncated 
the subtrees at their first node as they always have an odd number of leaves thereafter. 
The top left and bottom right nodes along the Mobius strip in (b) or (c) may also 
enter the lead for reflection quantities. 

this case) as the base diagram or structure. 

It is also beneficial to consult the full expansion in powers of the inverse channel 
number of the first two moments of the transmission eigenvalues [22, 25, 11, 26] and 
to draw the corresponding diagrams as ribbon graphs. The procedure of going from 
the closed periodic orbits to scattering trajectories and then to a graph is illustrated in 
figure 3 for the first subleading order correction. Removing the remaining subtrees from 
figure 3(f) leads to the base structure in figure 4(a) to which we can append subtrees 
to create valid diagrams like figure 4(b) whose boundary walk is depicted in figure 4(c). 
As the base structure involves a loop which is traversed in opposite directions by the 
trajectory and its partner, all the diagrams created in this way can only exist in systems 
with time reversal symmetry (corresponding to the orthogonal RMT ensemble). 

Along the loop we can add subtrees at any point and to make a valid /-encounter 
we must add 2/ — 2 subtrees (the remaining two stretches in the encounter belong to the 
loop itself). If the node has an odd number of trees both inside and outside the loop, 
we refer to it as an odd node. It is easy to convince oneself that in order to have each 
stretch of the loop traversed once by a 7-trajectory and once by a 7'-trajectory there 
must be an odd number of odd nodes around the loop. 

We start by evaluating the contribution of a node along the loop. For the node we 
include all possible sizes / of the resulting encounter. Adding the 21 — 2 subtrees (of 
which / — 1 start with an incoming direction and / — 1 with an outgoing direction) there 
are 2/ — 1 ways of splitting them into groups inside or outside of the encounter. With 
the / — 1 ways which result in an odd node we include the factor p whose power will 
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later count the total number of odd nodes around the loop. This leads to 

oo 

A{p) = j2Mffy-'[pii-i)+i]- (18) 

1=2 

The number of incoming and outgoing channels connected to the same node is equal 
[see an example in figure 4(b)]. An /-encounter can touch the lead only if every other 
edge connected to it is empty (connected directly to a leaf). Since for nodes on the loop 
we need to include the edges that belong to the loop itself and which cannot be empty, 
we conclude that only odd nodes can possibly touch the lead. Since in this case we need 
/ empty edges and we have / — 1 edges of each type, touching the lead is possible only if 
the incoming and outgoing channels are in the same lead, as they are when we consider 
a reflection quantity. 

With 2k — 1 trees on the inside of which k must be empty (and the remaining k — 1 
arbitrary) and the remaining 2/ — 2/c — 1 on the outside (with / — k empty and I — k — 1 
arbitrary) and with Zi = Zq = z for reflection quantities we add the following to the 
node contribution 

oo l—l 
1=2 k=l 

We then allow any number of nodes along the loop, though each time we add a new 
node it creates a new edge of the loop. Because of the rotational symmetry, we divide 
by the number of nodes. In addition, there is a symmetry between the inside and the 
outside of the loop, leading to a factor of 1/2. The total contribution thus becomes 

k=l 

Finally, to ensure that we have an odd number of odd nodes along the loop, we set 

This function then generates all the diagrams with 2n channels. We can now choose 
any of the leaves to be labelled ii, which fixes the numbering of all other leaves: they 
are numbered in order along the boundary walk. The freedom of choosing one of the 
leaves gives a factor of 2n. To get this factor we differentiate the result with respect to 
r and multiply by r so that the power of r still counts the total number of channels. 
Thus we obtain the generating function 

F = r^. (22) 
dr 

For the transmission, using the semiclassical values of the variables in (13), we find 
that the node contribution in (18) becomes 
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As we do not allow the nodes to enter the leads (as the incoming and outgoing channels 
are now in different leads) we also have B = 0. Note that the node contribution is given 
solely in terms of h = (fxp and the full contribution evaluates to 

^i = 7lnf^V (24) 



4 \l + h^ 

Putting in the correct explicit solution for h from (16) and transforming according to 
(22), we find the following generating function for the orthogonal correction to the 
moments of the transmission eigenvalues 

(1 - s)(l - s + 44s) 

where we set s = to generate the moments as the n-th moment involves 2n channels. 
This order correction was previously treated using a RMT diagrammatic expansion [9], 
and can be derived by performing an asymptotic expansion in inverse channel number 
of the RMT result for arbitrary channel number of [15]. 

For the reflection we have / = / and the node contributions in (18) and (19) are 

^(^) = — /i_/2)2 — ' B(p) = (l-rfy (^^^ 

Using relation (10) we can rewrite B{p) as 

pN p 

Ci (1 - / 

so that for the generating function we find 

where for the last term we simplified the numerator and denominator inside the 
logarithm by only keeping the remainder after polynomial division with respect to the 
quadratic for / in (11). Putting in the explicit solution from (11) and following (22) 
we obtain the rather simple generating function for the orthogonal correction to the 
moments of the reflection eigenvalues 

Note that this result only depends on ,^ = Ci(l — Ci) = CiC2; which is not so obvious from 
(28) and (11). However, the relations in (2) and the fact that the trace of the identity 
matrix, being the respective number of channels, is only leading order in inverse channel 
number means that the dependence only on ^ of the subleading transmission moments 
(25) must be mirrored in the reflection moments. For the reflection into lead 2 we simply 
swap and (2, which clearly does not affect this order correction. 
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(c) 1 (d) 




Figure 5. The first sublcading order semiclassical diagrams for systems witliout 
time reversal symmetry. We start witli tlie trajectory pairs wliicli contribute to tfie 
conductance in (a) and (b). Removing tlie clrannels and tlicir links we obtain the base 
structures (c) and (d) for this case. 

3.2. Unitary correction 

We can continue using the ideas above to treat higher order corrections. In particular, 
for systems without time reversal symmetry the first correction occurs at the second 
subleading order in inverse channel number. The semiclassical diagrams for the 
conductance are given, for example, in [26] and can be represented as the graph diagrams 
shown in figures 5(a) and (b). We note that this representation is not unique and 
is chosen for simplicity. It is also important to observe that, despite the twists, the 
corresponding ribbon graphs are orientable, i.e. have two surfaces (unlike the Mobius 
strip). It can be shown this is true in general: diagrams contributing to the unitary 
case are orientable. Further, the diagrams contributing at this order have genus 1, i.e. 
embeddable on a torus (but not a sphere). This, too, can be shown to continue: the 
contribution to the order l/A^^^^"^-* comes from diagrams of genus g. 

From the diagrams in figures 5(a) and (b) we can form the base structures by 
removing the channels and their links, see figures 5(c) and (d). A similar restriction to 
the one above still holds when appending subtrees to ensure that the resulting diagrams 
are permissible. Namely, the total number of odd nodes and twists along every closed 
cycle in the diagram has to be even. We note that the definition of an odd node depends 
on the cycle: the left node of figure 5(a) is odd with respect to the cycles formed from 
the top and bottom arcs and is even relative to the cycle formed from the top arc and 
the middle edge. We remark that this rule was enforced for the Mobius diagram as well. 

Finally, we need to discuss the symmetries of each base diagram. The generators 
of the symmetry group of base diagram 5c are the shift of the edge numbering and the 
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reflection, giving a group of size 6. Tlie generators for base diagram 5d are inside-out side 
mappings of the two edges, giving a group of size 4. 

When we append subtrees along each edge of a diagram, whether we append / or 
/-type subtrees depends in a complicated way on the types of subtrees appended along 
the other edges. Therefore we restrict ourselves to the simpler situation where / = / and 
only treat reflection quantities (with a = 0). Along the links connecting the nodes in the 
base diagrams we can append subtrees as before, but because the rotational symmetry 
is now broken we no longer divide by the number of nodes. The edge contributions can 
therefore be written as 

oo 

Eip) = y + B)f = YZT^^^y (30) 

where A{p) and -B(p) are as in (18) and (19) but with the simplification f = f and 
a = 0. 

We also need to append subtrees to the nodes of the base structures and, finally, 
ensure that we have the correct number of objects (odd nodes and twists) around each 
closed cycle. To proceed, we number each of the regions around the nodes and label the 
closed cycles with greek letters as in figures 5(c) and (d). We start with figure 5(c) and 
use powers of pa, Pp and p^ to count the number of objects along the respective cycles. 
At the top node we can add subtrees in any region we hke as long as we add an odd 
number in total to ensure that the top node becomes a valid /-encounter. An Z-encounter 
involves 21 stretches and we have 3 stretches already from the base structure. If we place 
ki subtrees in each region i and use the power of q to count the total number of subtrees 
added, we can write the contribution of the top node as 

oo 

with / = (/ci + A;2 + /cs -l- 3)/2 and where the number 3 in the subscript refers to the fact 
that the node in the base diagram starts with 3 stretches while the 'c' refers to its label 
in figure 5. Further, when we have an odd number of trees in each region, and when 
the odd numbered trees in each region are empty, then the top node can also enter the 
lead (since we are considering reflection quantities). If we define ki = 2ki + 1 then we 
have ki + 1 empty subtrees and ki arbitrary subtrees in each region. In total we would 
then add the contribution 



nM=QPc.PpP, (32) 

^1,^2,^3=0 



where I = {ki + k2 + k^ + 3) and we simplified the powers of q and p as in the end we 
are only interested in whether they are odd or even and they are all odd here. Finally 
to ensure that the total number of trees added is odd, we substitute 

V3,= ^^''^-^->'^--" +Vg,fa=l). (33) 
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The complete diagram in figure 5(c) is made up of two sucli nodes as well as three 
links. Each of the links lies on two cycles so we can write the full contribution as 

kl = ^-E{p^pp)E{ppp,)E{p,p^){y:,,f (34) 

where we divide by 6 to account for the symmetry of the base structure. Here the 
number in the subscript now refers to the order of the contribution while the 'U' in the 
superscript refers to the fact that these diagrams correspond to the unitary ensemble. 
Then to ensure that the number of objects along each cycle is even we simply average 

u kl{p = i) + kl{p = -i) 

^2c = 2 ' ^ ' 

for Pa, Pi3 and p^ in turn. 

For the base structure in figure 5(d) we now have a single node and four regions. 
Region 3 lies inside both cycles and each cycle starts with a single object inside (and 
likewise outside) which are the stretches of the other cycle leaving or entering the node. 
Treating the node as above, we obtain the contributions 

oo 

V,a{q)= MQff'^''^''^''¥a'{Pc.Pp)''p';\ (36) 

with I = (kl + k2 + h + k4 + 4) /2 and 

oo 

VM= E zif'^+'^+''+'^\ (37) 

where / = ki + k2 + k^ + k4^ + A and with an odd number of trees in each region in 
this second case we are guaranteed to add an even number to each cycle and an even 
number overall. With four edges touching the node in the base structure we need to 
add an even number of subtrees in total to the node to make a valid /-encounter so the 
node contribution reduces to 

V«^M^i>±^i^i^ + K.(,-l)^ (38) 
Including the two edges we have a total contribution of 

K^d = \paPisEipa)Eipfs)VM, (39) 

where we divide by 4 to account for the symmetry of the diagram and the PaP/s accounts 
for the fact that the cycles each start with a single object (the original node, which is 
odd for both cycles). We likewise take the average 

u k^p = 1) + k^,{p = -1) 

^2d = 2 ' ^ 

for Pa and p/j in turn. 

When we put in the semiclassical substitutions from (8) the formulae above can be 
summed and simplified. After applying the operator we find that first correction 
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Figure 6. The additional 4 base structures that exist for systems with time reversal 
symmetry at the second subleading order in inverse channel number. 



for the reflection for the unitary case (adding the two base cases) has the generating 
function 



By restricting ourselves above to the situation where / = / we are not able to obtain 
the transmission directly, but we can instead find the likely transmission generating 
function using (2) that tH + r\ri = I 



The fact that we get such simple functions is a little surprising especially as the result 
from each base case is notably more complex. In fact this pattern can be seen to continue 
if we expand the RMT result as in Appendix A. The generating function in (42) can 
also be obtained [15] from their RMT result. 

3.3. Second orthogonal correction 

When the system has time reversal symmetry, the edges and encounters can again 
be traversed in different directions by the trajectory set and their partners. For the 
conductance there are 7 further semiclassical diagrams at this order as depicted, for 
example in [26]. When we remove the starting and end links to arrive at the base 
structures, we find that they reduce to the 4 base cases depicted in figure 6. 

We note that the additional diagrams are non-orientable when viewed as ribbon 
graphs. Their groups of symmetry contain two elements each: reflection for diagrams 
6(a)-(c) and inside-out flipping of both edges simultaneously for diagram 6(d). 

Figures 6(c) and (d) are almost the same as figures 5(c) and (d), so we start 
with evaluating the contribution of figure 6(a). Although region 1 and 3 are spatially 




(41) 




(42) 
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connected they differ as to where we append subtrees at the nodes. Starting with the 
node on the left we therefore get the contributions 

oo 

V3M= E ^i{<lf)'^''^''^''^P^\ (43) 

ki,k2,k3=0 

with / = (/ci + + A^s + 3)/2 and 

oo 

VLiQ)=<lPo. J2 (44) 

fcl,fc2,fc3=0 

Again to ensure that an odd number of trees are appended, we substitute 

For the node on the right we obtain the contribution which is the same as Vsa but 
with Pa swapped with (and also /c2 with ki). 

Along with the two nodes in figure 6(a) we have three links, two of which form 
cycles which already contain a single object (a twist). The total contribution is then 

= \paP^E{pa)E{l)E{pp)V,S',^, (46) 
and to have an even number of objects along both cycles we average 

o kg{p = I) + kZ{p = -I) 

^2a = 2 ' *- 

for Pa and Pj3 in turn. 

The node in the base structure in figure 6(b) provides the following contributions 

oo 

VM = Yl x,(g/)('^+'=^+'^+'^V'^p'^ (48) 
with I = {ki + k2 + h + k4^ + 4) /2 and 

oo 

K{q)=Po.Pp zif'^+~''+~'^+~^^\ (49) 

The total number of trees added must be even, leading to 

T/ V^4b(g = i) + V^4b(g = -i) ^ , . . .... 

^4b = ^ + \/4b(g = l), (50) 

while with the two cycles (which each start with a single object) we have a total 
contribution of 

Kg = \paPpE{pa)E{pp)V,^. (51) 
As before we take the average 

o kg{p = I) + kg{p = -I) 

^2b = 2 ' ^ ' 

for Pa and P/3 in turn. 
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The difference of tlie structure in figure 6(c) from tliat in figure 5(c) is tliat now 
cycles a and /3 start witli an odd number of objects. Tlie contribution is tlien 

before averaging over tlie p's in turn. Only the a cycle in the structure in figure 6(d) 
now starts with an odd number of objects so its contribution is 

= ^p.E{p^)E{pp)V,^, (54) 

and then we average over p^ and p/3 in turn. 

Summing over the four new base cases (as well as the two that also exist without 
time reversal symmetry) we obtain the second subleading correction for the orthogonal 

case 

^^0(,) ^ _ ?.M3 + .) + l-2.1 ^ 
(l-4fs)i 

From this we can again find the likely generating function for the transmission 

NTOM ^ (56) 

(1 - 5)2(1 - s + 4^s)2 

The moments generated by (56) can be proven [15] to to agree with the moments 
obtained from an asymptotic expansion of their RMT result. Semi classically, time 
reversal symmetry allows more possible diagrams, so the results here are somewhat more 
complicated than the results (41) and (42) for systems without time reversal symmetry, 
but the RMT result [15] for the orthogonal case is notably more complex than for the 
unitary case. The results here are therefore useful in simplifying asymptotic expansions 
of RMT moments. 



3.4- Leading order revisited 

To obtain the leading order contributions [27], the start of the tree was fixed in the first 
incoming channel ii which allowed the top node to possibly i-touch. For example we 
simply place an incoming channel on top of the tree in figure 2(a). If the odd numbered 
subtrees after the top node were empty then the top node could z-touch the lead, leading 
to the generating function [28] 

00 00 

^0 = Cif + J2 + Coa + (7 5^ Zo,i{f + ay~\ (57) 
1=2 1=2 

Here the terms involving a derive from the fact that the section of the tree above the 
top node is actually an empty / tree. 

However, using the ideas we developed for the subleading corrections we can imagine 
a way of generating the leading order trees without fixing any of the channels as a root. 
As we shall see this is particularly beneficial for calculating energy dependent correlation 
functions as in sections 4 and 5. To start, we view a single point as the base structure for 
the leading order diagrams. We therefore obtain the leading order diagrams by joining 
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subtrees to this point. To create a valid encounter we need to add 21 subtrees (with 
/ > 2), / of which are /-type subtrees starting from an incoming direction and the other 
/ are /-type subtrees. If all of the subtrees of a particular type are empty then the node 
created can touch the lead. We obtain the generating function 

'^'""E^ + f^H-E^f^. (58) 

1=2 1=2 1=2 

where we divide by / because of the rotational symmetry and by a further factor of 2 
because of the additional possibility of swapping the incoming and outgoing channels. 
The last two terms in (58) represent moving the Z-encounter (formed by appending the 
subtrees to the starting point) into the incoming and outgoing channels. Importantly 
we overcount the trees by the factor V of their total number of encounters because, for 
a given resulting tree, any of the nodes could have been used as the base structure. 

On the other hand we can also construct trees by joining two subtrees together, one 
/-type and one /-type. After joining, the new vertex of degree two gets absorbed into 
the edge. A tree has exactly V — 1 internal edges, therefore there are V — 1 ways to obtain 
a given tree from joining two subtrees. The joining operation gives the contribution 

^0 = ^(/ - yCo){f - yCi) - ^CoiyCi + a), (59) 

where we divide by 2 because of the symmetry of swapping the incoming and outgoing 
channels. In the first term in (59) we subtract the empty tree from both / quantities to 
ensure that they both include at least one node so that the edge formed is an internal 
one. The last term in (59) is then to ensure that the diagram made of a single diagonal 
link with no encounters {V = 0) is included with the correct factor of —1/2. Taking 
the difference between (58) and (59) then means that we count each tree exactly once. 
Fittingly, for all the physical quantities we consider in this article (59) equals minus 
the / = 1 term in the sums in (58). This is natural because by joining two subtrees we 
essentially create a 1-encounter. We simplify the difference to 

Ko = Ko-K, = l^^^ + l^^ + l^ . (60) 

(=1 1=1 1=1 

Now that no root is fixed we can make any channel to be the first incoming channel 
and this generating function indeed misses a factor 2n compared to the generating 
function F. This turns out to be very useful for the density of states of Andreev 
billiards in section 4 and to recover F we can use relation (22). 

The generating function (57) for the reflection into lead 1 becomes 

^0 - (61) 

Taking the solution of (11) or inverting (61) and substituting into (11) leads to the 
generating function for the reflection 



Rojs) _ 2Cis - 1 + yi^^s 

N ~ 2(1 -s) • 
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For the reflection into lead 2 we swap (i and ^2- Note again that when we take the explicit 
solutions to any of the generating functions we chose the solution whose expansion in r 
agrees with the semiclassical diagrams. 

If, on the other hand, we start with (60), we obtain the generating function without 
the factor 2?7,, 

^ = ^ln(l-/')-Ciln(l-r/). (63) 

To obtain the missing factor 2n we substitute the correct solution of (11) into (63) and 
apply the operator r^, as in (22). Simplifying the result we recover (62). 

For the leading order transmission moments we also start with (61) so with (17) 
we obtain [27] 



iV 2 V 1 - s 2 ^ ^ 

The integrated generating function is 

^ = ^ ln(l - 00) - I ln(l - r0) - I ln(l - r0). (65) 

As for the reflection, substituting the solutions of (17) and transforming according to 
(22) we recover (64). 

4. Density of states of Andreev billiards 

If we imagine merging the scattering leads and replacing them by a superconductor, 
then our chaotic cavity becomes an Andreev billiard. Using the scattering approach 
[51], the density of states (normalised by its average) of such a billiard can be written 
as [52] 

(^^, (aa) 

n=l 

in terms of energy-dependent correlation functions of the full scattering matrix 



C{e,n) = j^TT 



sUe-^]s(e ■ 



(67) 



Here the energy difference is in units of the Thouless energy E'x = ^/2rd (which depends 
on the average classical dwell time Td) and measured relative to the Fermi energy E. 
Strictly speaking, for Andreev billiards we should use 5** (matrix S with complex- 
conjugated entries) instead of the adjoint matrix S'^ but we will only consider systems 
with time reversal symmetry where S is symmetric. With a superconductor at the lead, 
each time the particle (electron or hole) hits a channel it is retroreflected as the opposite 
particle (hole or electron) and semiclassically (see [29, 30] for fuller details) we traverse 
the partner trajectories in the opposite direction than for the reflection or transmission. 
For the leading order diagrams, this means that all the links (and encounters) are 
traversed in opposite directions by electrons and holes, so that if we break the time 
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reversal symmetry, with a magnetic field say, then none of these diagrams are possible 
any longer. Interestingly, at subleading order some diagrams are still allowed when the 
symmetry is completely broken, as for example the coherent backscattering contribution 
which comes from moving the node in figure 3(e) into the lead. 
We can consider the generating function 

oo 

G(s) = ^s"C(e,n), (68) 

n=l 

which generates the required correlation functions. We note that the definition of G here 
is marginally different than in [29, 30]. The semiclassical treatment there just requires 
us to make the substitutions 

y =J^^y x, = -N{l-la), 

Zi,i = Zo,i = r^N, Co = Ci = rN, a = 0, (69) 

where a = ie. 
4-1- Subtrees 

As Zj = Zo we have / = / so (5) becomes 

oo oo 

/(I - a) = r - ^(1 - la)^-' + r r''' f~\ (70) 

1=2 1=2 

or 

P2\ 



/(l-«-f 



(71) 



(l_/2)2 

which reduces to the cubic 

f-r{l + a)f-{l-a)f + r = 0. (72) 

4-2. Leading order 

The leading order in inverse channel number was obtained semiclassically in [29, 30], 
and for the energy dependent correlation functions we have 

^ f^. (73) 

1 - rf 

setting Go = F'q/N, inverting (73) and substituting into (72) we obtain the cubic 

(1 - sfGo^ + s{3s + a - 3)Go^ + s(3s + a - l)Go + = 0. (74) 
However, for the density of states the correlation function (60) 



f2 

N 2"^'^ ' ' ' 2(1 -~72^ 



^ = \ - f) + TTTl^ - ln(l - r/), (75) 



turns out to be more useful. Indeed with G = F/N and comparing (22) with (68), we 
see that 

— = y-C(6,n). (76) 

n=l 
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Introducing 

s"" dC{e, n) ^ d Ko 



n=l 

this is precisely what is required for the density of states of Andreev biUiards in (66), 
as by setting s = — 1 we have 

d{e) = l + 2ReH{s = -1). (78) 

Performing the energy differential implicitly, we arrive at the cubic 

a^{l - s)Ho^ + a [s{a - 2) + 2(1 - a)] Hq^ + [s{l - 2a) - (1 - af] Ho + s = 0. (79) 

Having the generating function K therefore allows us easier access to the density of 
states than we had previously [29, 30]. To make the connection to the RMT treatment, 
we set s = — 1 and make the final substitution 

SO that the leading order contribution to the density of states is do{e) = —lmWo{e) 
where Wq satisfies the cubic 

e^Wo^ + AeWo'^ + (4 + e'^)Wo + 4e = 0, (81) 

as found previously using RMT [53]. The result can be written explicitly as 

rfo(e) = Re^[Qo,+ (e)-Qo,4e)], (82) 

where 

Qo^^(e) = 8 - 36e2 ± 3eV3D '\ D = + Ue"^ - 16. (83) 



Note that the density of states is only non-zero when the discriminant D is positive, 
which occurs when e > 2 ^ • 

4-3. First correction 

With the techniques in this article we can go beyond this, and RMT, and see what 
happens at the next two orders in inverse channel number. For the first subleading 
order the even node contribution in (18) becomes 

P{ P-p-2) aP(P-3p + 2p + 4) 

:i - py {l-P 

while the odd node contribution in (19) is 

which we rewrite using (71). This simplifies the generating function (21), which reduces 
to 

- i [{l + a)f--i2 + a^)P + l-a) ' ^^'^ 
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We provide the further generating functions in Appendix B, but for the density of 
states we substitute Hi{s = —1) = iWi{e)/2 in (B.2) and obtain the cubic 

= -D\NWif + 2e(lle2 - 8)D{NWiy 

+ (e^ + 528e^ - 7206^ + 256)NWi + I6e{3e^ + 16), (87) 
so that the first correction to the density of states is given by 

di{e) = [Qi,+ (6) - gi,_(6)] , (88) 

with 

Qi^^(e) = [4096 - 46848e^ + 103584e^ - 3232e^ + 126e^ 

1 

± 3e(768 - 2928e2 + 96e^ + e^)V3D ' . (89) 
As this involves the same discriminant as (82), the correction to the density of states is 

a.o o„„ non.™ when . > . (^)'. .e. i. .Ke .a„. ,ap a. .Ke leadin, o.de. 
term. The correction however is negative and has a singular peak from the discriminant 
in the denominator. 

4.4- Second correction 

Repeating this procedure for the six base cases that contribute at the second subleading 
order we find 

d2{e) = -^Re^ [Q,A^) - Q^A^)] , (90) 

with 

g2,±(e) = [2eD^ (16224878592 + 740963778566^ + 153714421760e^ + 86120095744e^ 
+ 28154556672e* + 6522754176e^° + 739116528ei2 + 12120680e^^ 
+ 122837e^6 - 5324e^^) ± 3D (3189506048 - 416037601286^ 

- 187618951 168e^ - 192686981 120e^ - 237943482368e^ 

- 108211019264e^° - 13928492544e^2 + 1338160896e^^ 

1 

+ 274655180e^^ - 6774219e^® + 52756e2°) VSD ' . (91) 

This correction is positive again and has a larger and steeper singular peak than 
the first order correction. To illustrate this we plot the leading order result, as well as 
the two corrections in figure 7(a) for = 25, while in figure 7(b) we sequentially add 
the corrections to the leading order result, again for A^ = 25. 

The fact that the hard gap remains derives from the discriminant D which is already 
present in / (at s = = — 1) from (72). We could then expect that the gap is robust 
against further higher order corrections, seeing as the expressions always involve /. 
Someway above the gap we see that the corrections (especially the second) make little 
difference but it is the region directly above the gap which is particularly interesting. 
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Figure 7. (a) The leading order density of states (dotted) along with the first 
(dashed) and second (solid) correction for iV = 25. (b) The leading order density 
of states (dotted) with the first (dashed) and then the second (solid) correction added 
for N = 25. 



The expansion in inverse channel numbers is poorly (if at all) convergent but if the 
pattern of alternating singular peaks continues its sum (or rather the exact result for 
finite channel number) could take any value. In particular the gap could widen. 

Though treating the density of states of Andreev billiards semiclassically basically 
just involves using different values for the variables in the graphical recursions, on the 
RMT side it remains to be seen how one could use recent advances, like the Selberg 
integral approach, to proceed beyond the leading order in inverse channel number. For 
the leading order [53] a diagrammatic expansion was performed, but a result for an 
arbitrary number of channels would be especially welcome. It would determine whether 
the gap indeed persists and would clarify what happens to the density of states just 
above the current gap. 



5. Moments of the Wigner delay times 



Another quantity related to the energy dependent correlation functions are the moments 
of the Wigner delay times. To obtain them we define the correlation function 



D(e,n) 



N 



Tr 



(92) 



where we subtract the identity matrix in the form of S'^S = I. The corresponding 
generating function is 



L(.) = ^."D(e,n). 



n=l 



The moments of the Wigner-Smith matrix [54, 55] 
1 dE 



are then given by [28] 

fUn = Tr [Q]" 



i"n! de" " 



-D(e,n] 



e=0 



(93) 

(94) 
(95) 
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whose generating function we will denote 

oo „ 

Mis) = J2 ^^n. (96) 

n=l 

The identity matrix in (92) follows by removing the e dependence of the scattering 
matrices. However, as the identity matrix only has diagonal elements we identify these 
elements as diagonal trajectory pairs that travel directly from incoming to outgoing 
channels. We considered these to be formed when we moved encounters into the outgoing 
channels (equivalently we could use the incoming channels instead) whenever we formed 
an empty / subtree. The empty subtree is included in the general / contribution and we 
included the contribution a to allow us to change the effective value of / in this situation. 
The empty / subtrees consist of a single link and an incoming channel, producing the 
contribution yci = r/(l — a). To mimic subtracting the identity matrix in (92) we simply 
take away the value of this empty subtree at zero energy (a = 0) by setting 

a = -r, (97) 

along with the remaining semiclassical values in (69). 

5.1. Subtrees 

Including a breaks the symmetry of / and /, so the subtree recursions (5) and (6) 
become 

fjl-a-ff) ^ r /(I - a - //) ^ r 

(!-//> (l+r2-r/)' (!-//> (1-^/)' 

and we find [28] that / satisfies the following cubic 

(1 + r^f - r(l + r^){l + a)f - (1 + r^)(l - a)/ + r = 0, (99) 

while / is related by 

f = {l + r')f. (100) 

5.2. Leading order 

For the energy dependent correlation functions we have 
rNf r^N 
° ~ 1^7/ ~ (1 +r2)(l -r/)" ^^^^^ 

Setting Lq = Fq/N, inverting (101) and substituting into (72) we obtain the cubic 

(1 + s)Lq^ + as{l + s)Lq^ + s{2as + a - l)Lo + as^ = 0. (102) 

For the n-th moment of the delay times we want the coefficient of the n-th power of a 
which we can extract by transforming s ^ s/a and then setting a = 0. This leads to a 
quadratic and the leading order moment generating function [56, 28] 



1 - s - Vl - 6s + 



2 



Mo{s) = (103) 
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Alternatively we can start with the generating function from (60) 



N 2 ' y ■ ^ ■ 2(1 - (l + r2)/2) 

- ^ ln(l - r/) - i In [(1 + r'){l - rf)] (104) 

from which we can recover (102) and hence (103) by differentiating with respect to r, 
multiplying by r and using the result for ^ from differentiating (99) implicitly. 

5.3. First orthogonal correction 

For the first orthogonal correction we evaluate the contributions of the even and odd 
nodes around the Mobius strip 

_ hih-p-2) ahjh^ -3h + 2p + A) 

and 

Bip) = = (106) 

(l+r2-r/)(l-r/) (1 - h)^ ' 

which we again rewrite using (98) and both contributions only depend on h = ff. 
Putting these contributions into the generating function (21) we again obtain 

as in section 4.3 but with different values for / and / as given in (99) and (100). 
Differentiating in line with (22) and differentiating (99) implicitly we arrive at the 
generating function, given as (B.3) in Appendix B, which generates the orthogonal 
correction to the correlation functions D{e,n). Finally by transforming s ^ s/a and 
setting a = we find the correction to the moments of the delay times to be 



, . , , / N 1 — 3s — Vl — 6s + 
2(1 — OS + s'^) 

5.4. Next corrections 

Since we only treated reflection quantities where / = / for the next order corrections 
we cannot obtain the corresponding generating functions of the moments of the delay 
times. Instead we can generate the energy dependent correlation functions C(e,n) by 
expanding the generating function G2{s) which can be found by treating the six base 
cases as in section 4.4. Expanding to finite order, we can then obtain the functions 
D{e, n) using the relation 

n / 



k=0 



D{e,n) = }_^{-ir'l\C{e,k), (109) 



which follows from the binomial expansion of (92) and where C(e,0) = 1 has no 
subleading order contribution. Doing this only for the two base structures which exist 
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Figure 8. (a) With two traces, the semiclassical trajectories separate into two 
closed cycles. When the two sets do not interact as in (b) we recreate terms from 
Tr[XjXi]"i^ <^Tr[x|X2]"2^, but when they do interact as in (c) further diagrams 
are possible as in figure 9. 



without time reversal symmetry and plugging the resultant D{e,n) into (95) we obtain 
the moments for the unitary case to low order. We find that the generating function 

N'M^is) = f (110) 

(1 — OS + S^) 2 

fits with these moments, while if we treat all six base structures, the generating function 
iV^M?W = ^(i^H-5f(i^ill±^. (Ill) 

fits the low moments for systems with time reversal symmetry. 
6. Cross correlation of transport moments 

Along with transport moments, we can also consider non-linear statistics such as the 
cross correlation between transport moments, generated by 

oo 

Pix,M^i,s2)= sTsT (nMximxix.r) , (112) 

which involves two traces inside the energy average. Semiclassically we then have an 
expression involving two trajectory sets that form two separate cycles, as in figure 8(a). 
Of course when we look for trajectory sets which lead to a small action difference we can 
have independent diagrams for each set, as in figure 8(b). However, these are included 
in the individual moments treated previously, and when we remove them, 

00 

Plx.M'u S2) = P[xi,x.](5i, S2) - Yl ^r^'^ (rriXlX^r) (rriXlX^r) , (113) 

ni,n2=l 

we are left with trajectories that must interact, as in figure 8(c). In RMT, this quantity 
is known as the 'connected' part of the correlation function. 

It is interesting to consider the combinatorial interpretation of the interacting 
sets of trajectories. Denote the incoming channels belonging to the first trace by ij, 
j = 1, ... ,77-1 and the incoming channels from the second trace by i'j, j = 1, . . . ,n2- 
Input channels are mapped onto output channels by trajectories from X and also by 
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trajectories from . If we apply the first mapping followed by the inverse of the second 
mapping, we end up with the following transitions 

h-). Z2, ... z„jt->zi, and ■■■ ^ ^i- (H^) 

Thus the overall result is a permutation, written in the cycle notation as vr = 
{ii ... ini)ii'i ■■■ i'ni)- '^^^ interpret each Z-encounter as a cycle permuting / 
labels of the corresponding X trajectories. Then an entire leading order diagram can 
be interpreted as a factorization of vr into smaller cycles [c.f. [57]). In the language of 
combinatorics, interacting trajectories correspond to a transitive factorization, leading 
order corresponds to the minimality condition and the fact that encounters on different 
'branches' have no ordering imposed on them corresponds to counting inequivalent 
factorizations (up to a permutation of commuting factors). To summarize, the leading 
order interacting diagrams are in one-to-one correspondence with the minimal transitive 
inequivalent factorizations of a permutation into smaller cycles. This question has been 
studied combinatorially (for factorizations into transpositions only) for a permutation 
consisting of two cycles in [58] and for three and more cycles (these correspond to 3-point 
and higher cross correlations) in [59, 60]. We note that the above combinatorial questions 
and the evaluation of transport properties are not completely equivalent problems. 
To evaluate transport properties we need to find additional information regarding the 
number of encounters touching the lead. On the other hand, we make substitutions (8) 
or (13) which significantly simplify the results. 

6.1. Leading order 

We now proceed to expand the contributions of the interacting sets of trajectories in 
inverse powers of the total channel number by describing the corresponding graphical 
representations. The base diagram for the leading order term is just a single loop, like in 
figure 4(d) except with no twist. Without the twist, the two cycles of the permutation 
arise from the two walks on the inside and outside of the loop, see figure 9(a). One 
requirement, to ensure a small action difference, is that the parts of the loop are traversed 
on either side by parts of trajectories that contribute actions with different signs in the 
semiclassical expression. In figure 9 this means that the (blue) solid or dashed dotted 
lines on either side of the loop must partner (red) dashed or dotted lines on the other. 
Without time reversal symmetry the parts of the loop must additionally be traversed 
in the same direction by trajectory stretches and their partners, so that at odd nodes 
(those with an odd number of subtrees on each side of the loop) there is unequal number 
of channels of a given type, as in figure 9(a). With time reversal symmetry, parts of 
the loop may be traversed in any direction and we may also swap all the incoming and 
outgoing directions on one side (the inside say) of the loop as in figure 9(b). 

With these restrictions we can start to append subtrees at nodes around the loop. 
We will use the tree function fi and generating variable pi for the subtrees outside the 
loop and /2 and p2 for those inside. At each node we can either add an even or odd 
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Figure 9. Example graphs, which contribute to the leading order in inverse channel 
number of P[Xi.X2]i-^i^ ^2), made out of a single loop with subtrees attached at nodes. 
For a small action difference, the (blue) solid or dashed dotted lines must partner 
(red) dashed or dotted lines on the other side of the loop, (a) Without time reversal 
symmetry the odd nodes have either incoming or outgoing channels on the outer 
subtrees on each side of the loop, (b) With time reversal symmetry we may also 
swap the incoming and outgoing channels inside the loop. 



number of subtrees on each side of the loop and we start with the contribution when 
we add an even number 

1=2 k=0 1=2 ^ ^ 

where hi = /i/i. 

For systems without time reversal symmetry, with an odd number of subtrees on 
each side (an odd node) we have two possibilities as in figure 9(a). The subtrees in the 
odd positions on both sides all connect (first and last) to either incoming or outgoing 
channels. With outgoing channels we have the contribution 

<o = j:x,j2m-'f2ht'-'' =Y.^ihh \ , (116) 

1=2 k=l 1=2 ^ ^ 

while with incoming channels we swap / with / 

1=2 ^ ^ 

As before it is possible that an odd node touches the lead when the odd positioned 
subtrees on both sides are empty. Of course this also requires that the incoming or 
outgoing channels of the two quantities X\ and X2 originate or end in the same lead. 
When this is the case, we also have the following contribution 

00 i—\ 



r 

1=2 k=l 



00 



E 



■Plp2- 



1=2 P^f^ - P^f^ 
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where we needed to include the exphcit pi and p2 dependence of Zo on the number of 
empty trees on each side of the loop. With incoming channels instead we again swap / 
and / 

We allow an arbitrary number, /c, of nodes along the loop, but the total number of 
odd nodes must be even. Taking into account rotational symmetry, we get 

Aiia = -ln[l-yZ,(p)], (120) 

where Za(p) is a contribution of one node. The node can either be even or odd of one 
of two types. Since the number of odd z-nodes is equal to the number of odd o-nodes, 
we can write Z^{p) as 



ZM = ^'^^^ + P^iAlf + B^,){Aof + (121) 
To ensure that we indeed have an even number of odd nodes, we set 

leading to 

'^la = -I In [(1 - )2 - y2(^odd ^ 5^^^)(Aodd ^ 5^^^)] _ (123) 

For the additional contribution, in the case of systems with time reversal symmetry, 
from the diagrams like figure 9(b) we swap the incoming and outgoing channels inside 
the loop so that we now have the contributions 

<^ = T.^ifif2 \ _7 , <f = j:^ifih \ _7 . (124) 



hi-h2 ^ h-h 

and 



00 



^b,o - ^ -V— -T ' (125) 



1=2 ' Pl/l - P2/2 



00 



and correspondingly 

'^ib = In [(1 - )2 - y2(^odd ^ 5^^^^)(Aodd ^ _ (;L26) 

There is also an additional freedom of placing the label ii on any leaf outside, giving 
a factor of 2ni. Once ii has been placed the type (in- or out-) of the leaves inside is 
fixed and the freedom of placing the label i'l inside produces only a factor of n2. We 
obtain these factors by differentiating with respect to the variables pi and p2, 

r = (127) 

2 dpidp2 
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We further note that the differentiation ensures that there are at least two channels 
both inside and outside the loop. 

Using the appropriate semiclassical values of the variables £C, z and y as well as the 
corresponding subtree contributions, we find the following generating functions 

1 -2e(si + s2) 



-P[ri,ri],l('5l' 



-^[ri,r2],l('^l' -^2) 



S\S2 
2(Si-S2) 



- 1 



(128) 



with S\ = p\ and S2 = p\ and twice this result for the orthogonal case with time reversal 
symmetry. Even though for the autocorrelation P^^ ^ we can always move the odd 
nodes into the lead while for the cross correlation P^^ ^ we can not, this is somehow 
compensated for by the different subtree contributions and both give the same result. 

For the transmission autocorrelation we can move the odd nodes into the lead only 
for the unitary diagrams, leading to the generating function 



S1S2 



2(S1-S2) 



1 + (2e - 1) (si + S2) 



40 S1S2 



1-^2^/1 



(129) 



and we still obtain twice this for the orthogonal result. Finally for the cross correlation 
between the reflection and transmission we have 



-Pln.tl.ll-^l' -52) 



2 (si - S1S2 + S2) 



1 - S2 + 2^ (S2 - Si + S1S2) 



v/i-4esi(i-s2)ji 



1-S2 J 



(130) 



and twice this for the orthogonal case. Note that the above results remain unchanged if 
we swap ri and r2 as this just means swapping C,\ and C2 in the semiclassical contributions 
which does not change ^. From these results we can obtain the corresponding P\Xx,x-i\ up 
to the first subleading order by including the first three orders in inverse channel number 
of the moments corresponding to X\ multiplied by the moments corresponding to X^. If 
we also include n = terms (which are just the number of channels in the respective lead) 
with those moments then we obtain the ni = and ^2 = terms in (112). This then 
allows us to check that expansions of the various transport correlation functions indeed 
fulfil the unitarity conditions in (2). Note that if we assume a 'priori that the unitarity 
is preserved by the semiclassical approximation (3), any one of equations (128)-(130) 
implies all others. 



6.2. Subleading correction 

We can continue this process and look at the base structures like in figures 5 and 6 
but which separate into two cycles. In fact the possibilities are almost the same as in 
figure 6 but with one twist more or fewer as depicted in figure 10. These can also exist 
only for systems with time reversal symmetry and we can treat them in a similar way 
as before, but with the modifications above. 
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(d) 




Figure 10. The base structures which break into two cycles (for systems with time 
reversal symmetry) at the second subleading order in inverse channel number. 



Again the types of subtrees at each node depends on the nodes elsewhere in the 
diagram, so we restrict out attention to the simpler case of the reflection where / = /. 
Because of (2), we have 



Tr[r|ri]" - Ni 



Ti[thf = Ti[rlr2r - N2, 



(131) 



so that Pu 



P[ri,r2] where the reflection autocorrelation is equal to the reflection 
cross correlation. However for the cross correlation, as the channels of the two reflections 
are in different leads, the nodes which lie on both cycles can not enter the lead, and 
this further simplifies the calculation. The edges which travel through both cycles then 
provide factors 

1 



v{p) 



(132) 



1 - (A'^vcn _^ pylodd) ' 

while the edges which only pass through one cycle provide factors E[p) as before, but 
with the semiclassical values corresponding to the reflection into lead 1 or lead 2 as 
appropriate. We will denote this correspondence by a subscript in the following. 

The treatment of the diagrams is very similar to that in section 3.3 so we merely 
highlight the steps here. But first we discuss the symmetry factors. Because of time- 
reversal symmetry, we can put the first incoming channels on both faces on any leaf, 
leading to the differential operator pip2 Qp^Qp^ ■ Unlike the diagrams in figure 9, the 
'inside' and 'outside' faces are, in general, not related by symmetry, and we should 
consider both putting /i-trees on the 'outside' and on the 'inside'. For brevity we will 
only list the former contributions. Finally, the symmetry groups of diagrams 10(c) and 
10(d) have order 2 and 4 correspondingly and we will divide their contributions by the 
appropriate factor. 
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Starting with the diagram in figure 10(a), for the node on the left, which cannot 
enter the lead, we have 

oo 

^3a(g) = - i(lht'^''\<lhf'pl\ (133) 

fcl,fc2,fc3=0 

with an odd number of trees appended 

For the node on the right we have as before but with semiclassical values 
corresponding to the refiection into lead 1. To ensure a valid semiclassical diagram 
we still need each cycle to contain an even number of objects, so we have 

^^2! = PpVi.Pa)Ei{X)Ei{pp)v^J/^^. (135) 

This is then averaged 

,0, ^?.(P-1)+^?.(P--1) _ (13,) 

for Pq, and in turn. Finally, we add the contribution where we place /i subtrees along 
the inside and /2 subtrees along the outside. For the reflection cross correlation this 
reduces to swapping pi with p2 and swapping with ^2 = 1 — Ci- 
The node in figure 10(b) gives 

00 

^4b(g) = - i<iht'^''^''K<if2f'pi%'^ (137) 

^1,^2,^3,^4=0 



with an even number of subtrees in total 

2 



t>4b(g = 1) + VAh{q = -1) ,^ooN 



The total contribution is then 

'^"b = Pl3ViPa)Ei{p(3)l^ih- (139) 

averaged over pa and Pi3 in turn and we again add the result where we swap the trees 
on the inside and the outside. 

For the structure in figure 10(c) the nodes provide 

00 

hc{q) = - i<lhf'i<lht'^''^Pi%'p';\ (140) 

with 

so that the contribution is then 

k2c = \pc^Pp E2 {PaPp ) V {P^Pi ) V {PiPa ){v^cf (142) 

before averaging over the p's in turn. Likewise we include the contribution where we 
swap the subtrees on the inside with those on the outside. 
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Finally the node in figure 10(d) provides 

oo 

= - E (?/i)^'^"''^H9/2)('^+'=^V^^(p.P/.)'=^P^ (143) 

with 

= ^^'^(^ = ^) +^^^^(^ = _ (144) 
The total contribution is then 

l^2d = \v{PaMPl3)l'id, (145) 

averaged over and in turn. We also include the contribution where we swap the 
subtrees on the inside and outside. 

Summing the four diagrams, we eventually arrive at the generating function 

NP[?„r,],2i^u S2) = [si + S2- 2siS2 - 2^ (4s? + + 3siS2 -si- 5sls2 - 2s,sl) 

+ 8sl^^ {si + 3s2 - 3siS2 - si)] 

X 7 N3n T/ ^2 /T- ^r= + ^ "2) , (146) 

[Si - S2y{l - 4^si) VI - 4^S2 
where (si O S2) means we add the result with si and S2 swapped. 

We could check that the results in this section all agree with the first four moments 
calculated from the arbitrary channel RMT results of [12]. 



7. Conclusions and discussion 



We described a method for the semiclassical calculation of the expansion of several 
transport statistics asymptotically in the inverse channel number 1/A^. The calculation 
is performed by grafting trees onto the base structures with a low number of cycles, 
and relies on the fact that attaching trees does not change the order in inverse channel 
number. Instead, the trees add more incoming and outgoing channels and so increase 
the order of the moment. With graphical recursions, this allows us to generate all 
the moments at a given order in inverse channel number, which we performed up to 
the third order. The terms we considered suggest the following observations about the 
ribbon graphs that arise as the contributing diagrams 

• absence of time reversal symmetry results in graphs being orientable; both 
orientable and non-orientable graphs contribute to the calculation with time reversal 
symmetry, 

• the order (in 1/A^) of a contribution is reflected in the genus of the corresponding 
graph, 

• linear moments (with one trace) result in graphs with one face, while non-linear 
moments with m traces will require considering graphs with m faces. 
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The above general observations suggest that a complete expansion in should be 
both feasible and interesting to specialists working in algebra and combinatorics. 

We find that the semiclassical contribution of individual base diagrams depends 
significantly on the global structure of the diagram. This is in contrast to the expansions 
of the first two moments performed in [25, 11, 26], where the contribution factorized 
into a product over the vertices of the diagram and the problem was thus reduced to a 
combinatorial enumeration. The latter was achieved through finding recursion relations 
which connect different diagrams, and similar ideas could well be useful for the base 
structures we need for all moments. To illustrate the scale of the problem of going 
to higher order in 1/A^, for the unitary case there are 1848 base diagrams at the next 
contributing order. As they can involve more than two nodes, they can no longer be 
derived by cleaning the corresponding semiclassical conductance diagrams as was the 
case for the orders treated in this article. However, in the end all these diagrams would 
probably lead to the generating functions (A. 3) and (A. 7), highlighting the scale of the 
simplifications that take place. 

Our results fully agree with the predictions of RMT theory (as far as those are 
available [15]), and importantly are given in terms of very simple generating functions. 
This would suggest that extending the types of asymptotic analyses of [47, 48, 49] 
beyond the leading order, (as is currently being performed [15]), one could also expect 
to see simplifications of the RMT results. For the unitary case, where several different 
formulae are known for the moments of the transmission eigenvalues [13, 14, 15, 16], 
this analysis and the semiclassical endpoints could shed light on the combinatorial 
relationships between the different approaches. Because of the connection between 
RMT and weakly disordered systems, we can expect that our results also apply to such 
systems. Likewise, with the close correspondence between semiclassical and disorder 
diagrams [43], one might hope to find similar graphical recursions in a perturbative 
expansion of the appropriate nonlinear a model. 

If we were to consider non-linear statistics involving three traces (with their 
mean parts removed), the leading contribution would come from the diagrams like in 
figures 6 and 10 which split into three cycles, i.e. diagrams (a)-(c) without any twists. 
Considering quantities with m traces we would then need to treat base diagrams related 
to those which contribute to order l/A^(™-~^) (and higher) for the linear statistics, and so 
we immediately run into the considerations and difficulties described above. Curiously 
though the moments of the conductance and the shot noise themselves can be efficiently 
treated using RMT [14, 17, 18, 19]. Semiclassically, the m-th moment of these quantities 
corresponds to having exactly 2 or 4 channels along each of the m cycles, and the RMT 
results might then provide a pathway for generating such semiclassical diagrams. This 
could in turn be useful for generating and treating the corresponding base diagrams for 
the linear statistics. 

The methods described in this article were also used to treat the density of states 
of Andreev billiards. Replacing the normal conducting leads of the chaotic cavity by a 
superconductor produces strong effects like the complete suppression of the density of 
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states around the Fermi energy [53, 29, 30]. Being interested in the density of states 
one must evaluate moments of all orders, something that our methods are particularly 
geared towards. Going beyond leading order in inverse channel number we could show 
that this gap persists for the next two orders, and that the behaviour of the density of 
states slightly above the gap is not determined by just these terms in the expansion. 
Because of the superconductor, one not only needs to know all the moments but also 
all the higher orders in inverse channel number. A result for arbitrary channel number 
would therefore be particularly welcome for such systems, which leads to the question 
of how to adapt the recent RMT advances to tackle this problem. Similarly, chaotic 
cavities with additional superconductors attached (Andreev dots) also exhibit significant 
effects due to the presence of the superconductors and also require one to be able to 
treat what would correspond to all the moments of usual transport quantities. For 
example, at leading order in inverse channel number, the conductance through a normal 
chaotic cavity requires just the diagonal pair of trajectories, while for the conductance 
through an Andreev dot one needs full tree recursions [32]. The treatment is actually 
similar to the edges in the base diagrams here, but with the added ingredient of having 
two (or more) different species of subtree. One can then see that treating the transport 
moments of Andreev dots requires an extra layer of complexity compared to normal 
chaotic cavities. 

The results in this article are all for the case in which the leads are perfectly coupled 
to the chaotic cavity, rather than for the more general and experimentally relevant case 
of non-ideal coupling. This is typically modelled by introducing tunnel barriers into 
the leads with some probability to backscatter when entering (or leaving) the cavity. 
Semiclassically, along with affecting the contributions of the channels and modifying the 
survival probability and hence contributions of the links and the correlated trajectory 
stretches inside the encounters, the main change is that a wealth of new diagrams become 
possible [61]. Specifically, encounters may now partially touch the leads and have some 
of their links backreflected at the tunnel barrier while the rest tunnel through to enter 
or exit the system. In principle, these possibilities would become extra terms in the tree 
and graphical recursions in this article, but so far these types of diagrams have only been 
treated semiclassically for the lowest moments [61, 62]. However, from a RMT viewpoint 
at leading order in inverse channel number and with the same tunnelling probability for 
each channel, the non-ideal contacts just increase the order of the generating functions 
by one, for example for both the density of states of Andreev billiards [53] and the 
moments of the Wigner delay times [63]. 

Finally one can wonder whether the effect of the Ehrenfest time can be incorporated 
into the graphical recursions developed here. For the leading order in inverse channel 
number, the effect could be included [33] in the tree recursions. First, the trees are 
related to each other through a continuous deformation, for example by giving the nodes 
a certain size (actually, the Ehrenfest time itself) and allowing them to slide into each 
other. Second, this is then partitioned in a particular way so that one can extract the 
Ehrenfest time dependence efficiently. Each partition and hence the sum of all diagrams 
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leads to the same simple Ehrenfest time dependence at leading order in inverse channel 
number [33] and this pattern and treatment seems to also hold at the first subleading 
order [64]. Whether this continues to higher order and nonlinear statistics, which start 
to include the complications of periodic orbit encounters, is an intriguing question. 
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Appendix A. Unitary reflection and transmission 



Using the RMT result from [14] we can compute the moments up to finite order and 
expand them in powers of the inverse channel number. Looking at the patterns for the 
reflection in (62), (29) and (41), we can expect that each order in the inverse channel 
number just increases the powers in the denominators (the square root comes simply 
from the subtrees) . In fact we find that the first several subleading orders can be written 
as 
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Similarly we can write the transmission as 
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V 

Similar patterns hold for the moments of the delay times for the unitary case, and 
the results are actually simpler than for the transmission and reflection since there is 
one parameter fewer. The likely generating functions can be found by expanding the 
RMT result of [15] and fitting to the behaviour of (103) and (110). 



Appendix B. Further generating functions 



For the energy dependent correlation functions, we find that the generating function 
Gi = Fi/N satisfies the cubic 

= [4(1 - af + s(a^ - 200^ - 8) + 4s{a + if] ^ {NGif 

+ 2 [{a- if + s2(a + 1)^] [4(1 - af + s{a^ - 20a^ - 8) + As'^{a + if] {NGif 
+ [{a - if - 4s(a + l)(a - 1)^ + 3s'^{a - l)(a + l)(a^ - 8a^ - 2) 

+4s^{a - l)(a + 1)^ + s\a + 1)^] NGi 
+ a'^s[{a-lf + s^{a + lf], (B.l) 

while the energy differentiated generating function instead satisfies 
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X [4(1 - + s{a^ - 20a^ - 8) + As^{a + if] (iVifi)^ 
+ [{a - 4)2(a - 1)^ + 4s(a - lf{6a^ - 21a^ + 4a - 16) 
+ s^{-3a^ + 156a^ + 102a^ + 96) 

-4s^(a + l)2(6a^ - 210^ + 4a - 16) + /(a + Af{a + if] NHi 

+ as{s - 1) [{a + 2)(a - + s(a - 2)(a + 4)^] . (B.2) 

For the moments of the delay times we have the first subleading order generating 
function 

= (s + 1)2 [4(1 - a)3 + sa{a^ - 8a^ + 4a - 24) + s'^a^{a^ + 4)] ^ (iVLi)^ 
+ 2(s + 1) [(a - 1)^ + 2s(a - 1)^ + 2s^a{a^ + 3)] 

X [4(1 - af + sa{a^ - 8a^ + 4a - 24) + s'^a^{a^ + 4)] {NLif 
+ [(a - 1)^ + 4sa(a - 3)(a - 1)^ + 3s'^a'^{a - l)(3a^ - 13a2 + 20a - 28) 

+2s^a2(a - l)(5a^ - lla^ + 16a - 32)^ + s^a2(5a^ + 27a2 + 32)] NLi 
+ sa^ [(a - 1)^ + 2s(a - 1)^ + 2s^a{a^ + 3)] . (B.3) 
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